home *** CD-ROM | disk | FTP | other *** search
/ POINT Software Programming / PPROG1.ISO / pascal / swag / math.swg / 0081_Getting Big Pi!.pas < prev    next >
Encoding:
Pascal/Delphi Source File  |  1995-03-26  |  11.9 KB  |  455 lines
  1. (*                          PROGRAM PI5.PAS
  2.                            By  David Adamson
  3.                            October  23, 1994
  4.  
  5. This it the third version of a program for computing Pi up to 2150
  6. places.
  7.  
  8. The first version was written in Turbo Pascal 3. The second was also
  9. in Turbo Pascal, but with all the math procedures converted to inline
  10. for maximum program executing speed.
  11.  
  12. This version is rewritten to allow compilation with Turbo Pascal 5 and
  13. up.
  14.  
  15. The inline code was scrapped for version 3. The sacrifice in speed is
  16. not a problem, because it improvement in today's computers overshadows
  17. the processing efficiency loss.
  18.  
  19. Here are my results using a 33MHz 386.
  20.  
  21. Pi v2                               Pi v3
  22. 100  places =  0.1 sec              100  places =   0.3 sec
  23. 1000 places = 11   sec              1000 places =  32   sec
  24. 2135 places = 48   sec              2135 places = 150   sec
  25.  
  26.  
  27. ----------------------------------------------------------------------
  28. THIS IS THE ORIGINAL HEADER:
  29.  
  30.                             PROGRAM PI1.PAS
  31.                               TANDY 2000
  32.                            BY CINO  HILLIARD
  33.                             APRIL 6 , 1986
  34.        ( AN IMPROVED VERSION OF MY ORIGIONAL PI.PAS IN THE  DL)
  35.  
  36. THIS PROGRAM COMPUTES THE DIGITS OF PI USING THE ARCTANGENT FORMULA
  37.  
  38.                 (  TERM 1  )   (  TERM 2  )
  39.     (1)  PI/4 = 4 ARCTAN 1/5 - ARCTAN 1/239
  40.  
  41.                                    IN CONJUNCTION WITH THE GREGORY SERIES
  42.                      N
  43.     (2)  ARCTAN X = SUM [ (-1)^M*(2M + 1)^-1*X^(2M+1) ]  (APPROXIMATELY)
  44.                     M=0
  45.  
  46. SUBSTITUTING INTO (1) AND (2) THE FIRST FEW VALUES OF N, REARRANGING,
  47. SIMPLYIFING, AND NESTING  WE HAVE,
  48.  
  49. PI= 3.2 + 1/25[-3.2/3 + 1/25[3.2/5 + 1/25[-3.2/7 + ...].].]       ( TERM 1 )
  50.     -1/239[4 +1/239^2[-4/3 +1/239^2[4/5 +1/239^2[-4/7 +...].].]   ( TERM 2 )
  51.  
  52. USING THE LONG DIVISION ALGORITHM AND SOME TRICKS I DISCOVERED, THIS (
  53. NESTED ) INFINITE SERIES CAN BE USED TO CALCULATE PI TO A LARGE NUMBER
  54. OF DECIMAL PLACES IN A REASONABLE AMOUNT OF TIME.  A TIME FUNCTION IS
  55. INCLUDED TO SHOW HOW SLOW THINGS GET WHEN N IS LARGE.
  56.  
  57. IMPROVEMENTS CAN BE MADE BY INCREASING THE NUMBER OF DIGITS EACH ARRAY
  58. ELEMENT HOLDS AND BY CHANGING THE DATA TYPE FROM INTEGER TO REAL FOR
  59. SELECTED VARIABLES. OF COURSE THE ADDED NUMBER OF DIGITS THESE CHANGES
  60. PRODUCE WILL COST MUCH MUCH MORE TIME. AH INDEED, 'TIS NO FREE LUNCH!
  61. HOWEVER, SINCE TERM 1 AND TERM 2 ARE COMPUTED SEPERATELY AND SINCE THE
  62. ARRAYS ARE STEP BY STEP UPDATED, THE  PROGRAM DOES LEND ITSELF TO
  63. PARALLEL OR NON VON ( WHAT A COINCIDENCE - SEE BELOW )  PROCESSING.
  64.  
  65. FOR EXAMPLE, LET COMPUTER 1 PERFORM TERM 1 AND COMPUTER 2 PERFORM
  66. TERM 2. MOREOVER, LET SEVERAL COMPUTERS SHARE IN THE CALCULATION OF
  67. EACH OF THE INDIVIDUAL TERMS. HOWEVER, TO KEEP EACH COMPUTER EQUALLY
  68. BUSY, A LOGARITHMIC TYPE OF ADJUSTMENT MUST BE MADE TO DECIDE ON THE
  69. NUMBER OF TERMS TO BE ASSIGNED TO EACH COMPUTER. SINCE THE HIGHER
  70. POWER TERMS REQUIRE MUCH LONGER TO COMPUTE, THE COMPUTERS ASSIGNED TO
  71. HIGHER POWERS MUST BE GIVEN  FEWER TERMS TO DO.
  72.  
  73. A LITTLE HISTORY
  74. ----------------
  75. IN AUGUST, 1949, PROFESSOR JOHN VON NEUMANN USED FORMULAS  (1) AND (2)
  76. TO CALCULATE PI TO 2035 DECIMAL PLACES ON THE  ENIAC  COMPUTER. THE
  77. EFFORT WAS MADE TO DETERMINE IF THE DIGITS CONFORMED TO SOME TYPE OF
  78. PATTERN OR IF THEY WERE RANDOMLY DISTRIBUTED. THE CALCULATION WAS
  79. COMPLETED OVER THE LABOR DAY WEEKEND WITH THE COMBINED EFFORTS OF FOUR
  80. ENIAC STAFF MEMBERS WORKING IN EIGHT-HOUR SHIFTS TO ENSURE CONTINUOUS
  81. OPERATION OF THE ENIAC. THE CALCULATION (INCLUDING CARD HANDLING TIME)
  82. TOOK APPROXIMATELY 70 HOUR.
  83.  
  84. THE CONCLUSION WAS AS SUSPECTED - THE DIGITS APPEARED TO BE RANDOM!
  85.  
  86. SOME YEARS AGO I REQUESTED INFORMATION ON PI FROM THE ENCYCLOPEDIA
  87. BRITANNICA RESEARCH SERVICE. I RECEIVED A REPORT GIVING THE ABOVE
  88. HISTORICAL ACCOUNT PLUS A LISTING OF THE 2035 DIGITS.
  89.  
  90. A LITTLE COMMENT
  91. ----------------
  92. USING MY T2K, PI1.PAS COMPUTES PI TO 2035 PLACES IN  15 MIN  31.97
  93. SEC.
  94.  
  95. GEE WHIZ!  WHAT EARTHLY USE CAN BE MADE OF THIS?  I PLAN TO MAKE
  96. DESIGNS AND COLOR PATTERNS USING THE DIGIT VALUES AS THE BASIS. MAYBE
  97. SOMETHING LIKE FRACTALS WOULD BE A GOOD START. THIS IS WHY COMPUTATION
  98. AND SPEED ARE IMPORTANT AS I DO NOT WISH (NOR TRUST) TO KEY IN ENTRIES
  99. FROM A PUBLISHED LIST.
  100.  
  101.                              CINO HILLIARD
  102.                               [72756,672]
  103.  
  104. ---------------------------------------------------------------------------
  105. MODIFICATIONS MADE FOR VERSION 2:
  106.  
  107. ADDED INLINE CODE TO IMPROVE THE SPEED. USING THE ORIGINAL SOURCE ON
  108. MY TANDY-1000  ( WITH A V20 CPU ) , COMPUTING PI TO 2035 PLACES RAN IN
  109. 39 MIN 44.92 SEC.. AFTER ADDING THE INLINE CODE IT RAN IN 10 MIN 25.16
  110. SEC..
  111.  
  112. ADDED A TIME FUNCTION FROM 'TECH JOURNAL FEB 85'. THIS USES A MSDOS
  113. INTERUPT, SO IT SHOULD RUN ON ANY MSDOS COMPATIBLE COMPUTER.
  114.  
  115.                               CHUCK WHITE
  116.                              [75006,3677]
  117. *)
  118.  
  119. PROGRAM PI5;
  120. Uses DOS, CRT;
  121.  
  122. Type
  123.   TimeString = string[12];
  124. VAR K,I,I2,J,M,N,Q,V,R,D,Z : INTEGER;
  125.     A,P,T : ARRAY[0..5000] OF integer;
  126.     TI,T2 : STRING[20];
  127.  
  128.  
  129. function Time: TimeString;
  130.  
  131. var
  132.    Hour, Minute, Second, Sec100 : word;
  133.    Hr, Min, Sec, Hun            : string[2];
  134.  
  135. begin
  136.    gettime(Hour, Minute, Second, Sec100);
  137.    str(Hour:2, Hr);
  138.    str(Minute:2, Min);
  139.    str(Second:2, Sec);
  140.    str(Sec100:2, Hun);
  141.    if Hr[1]  = ' ' then Hr[1]  := '0';
  142.    if Min[1] = ' ' then Min[1] := '0';
  143.    if Sec[1] = ' ' then Sec[1] := '0';
  144.    if Hun[1] = ' ' then Hun[1] := '0';
  145.    time := Hr+ ':'+ Min+ ':'+ Sec+ '.'+ Hun
  146. end;
  147.  
  148. PROCEDURE DIV32IA;       { DIVIDE 3.2 BY I AND STORE IN ARRAY A }
  149.  BEGIN
  150. Q:=3 DIV I;
  151. R:=3 MOD I;
  152. A[0]:=Q;
  153. V:=R*10+2;
  154. Q:=V DIV I;
  155. R:=V MOD I;
  156. A[1]:=Q;
  157. FOR J:=2 TO M DO
  158.       BEGIN
  159.       V:=R*10;
  160.       Q:=V DIV I;
  161.       R:=V MOD I;
  162.       A[J]:=Q;
  163.     END;
  164.  END;
  165.  
  166. PROCEDURE DIVA(D:INTEGER); { DIVIDE A BY SPECIFIED INTEGER FROM }
  167.  BEGIN                     { PROCEDURE COMPUTE; }
  168.     R:=0;
  169.     FOR J:=0 TO M DO
  170.      BEGIN
  171.      V:= R*10+A[J];
  172.      Q:= V DIV D;
  173.      R:= V MOD D;
  174.      A[J]:=Q;
  175.      END;
  176.  END;
  177.  
  178. PROCEDURE DIV4IA;         { DIVIDE 4 BY I AND STORE IN A }
  179.  BEGIN
  180. Q:=4 DIV I;
  181. R:=4 MOD I;
  182. A[0]:=Q;
  183. FOR J:=1 TO M DO
  184.      BEGIN
  185.      V:= R*10;
  186.      Q:= V DIV I;
  187.      R:= V MOD I;
  188.      A[J]:=Q;
  189.      END;
  190.  END;
  191.  
  192.  
  193. PROCEDURE DIV32IP;        { DIVIDE 3.2 BY I AND STORE IN P }
  194.  BEGIN
  195. Q:=3 DIV I;
  196. R:=3 MOD I;
  197. P[0]:=Q;
  198. V:=R*10+2;
  199. Q:=V DIV I;
  200. R:=V MOD I;
  201. P[1]:=Q;
  202.    FOR J:=2 TO M DO
  203.      BEGIN
  204.      V:= R*10;
  205.      Q:= V DIV I;
  206.      R:= V MOD I;
  207.      P[J]:=Q;
  208.      END;
  209.  END;
  210.  
  211.  
  212. PROCEDURE DIV4IP;       { DIVIDE 4 BY I AND STORE IN P }
  213.  BEGIN
  214. Q:=4 DIV I;
  215. R:=4 MOD I;
  216. P[0]:=Q;
  217.     FOR J:=1 TO M DO
  218.      BEGIN
  219.      V:= R*10;
  220.      Q:= V DIV I;
  221.      R:= V MOD I;
  222.      P[J]:=Q;
  223.      END;
  224.  END;
  225.  
  226. PROCEDURE SUBA;         { SUBTRACT A FROM P  AND STORE IN A }
  227.  BEGIN
  228.    FOR J:=0 TO M DO
  229.    A[J]:=P[J]-A[J];
  230. END;
  231.  
  232. PROCEDURE SUB32A;        { SUBTRACT A FROM 3.2 AND STORE IN T }
  233.  BEGIN
  234.    T[0]:=3;
  235.    T[1]:=1;
  236.    FOR J:=2 TO M DO
  237.    T[J]:=9-A[J];
  238.  END;
  239.  
  240. PROCEDURE SUB4A;        { SUBTRACT A FROM 4 AND STORE IN A }
  241.  BEGIN
  242.    A[0]:=3;
  243.    FOR J:=1 TO M DO
  244.     A[J]:=9-A[J];
  245.  END;
  246.  
  247. PROCEDURE SUBT;         { SUBTRACT TERM2 FROM TERM 1 AND STORE IN A }
  248.  BEGIN
  249.    FOR J:=M DOWNTO 1 DO
  250.    BEGIN
  251.     A[J]:=T[J]-A[J];
  252.     WHILE A[J]<0 DO
  253.      BEGIN
  254.       A[J]:=A[J]+10;      { ADJUST FOR NEGATIVE }
  255.       A[J-1]:=A[J-1]+1;
  256.      END;
  257.       WHILE A[J]>9 DO    { ADJUST A FOR EXCESS CARRY }
  258.        BEGIN
  259.         A[J]:=A[J]-10;
  260.         A[J-1]:=A[J-1]-1;
  261.        END;
  262.    END;
  263.  END;
  264.  
  265.  
  266. PROCEDURE COMPUTE;      { COMPUTE THE NESTED  SERIES  FOR I2 ITERATIONS }
  267. BEGIN
  268.  I:=I2;                  { COMPUTE TERM 1 TO I2 PLACES }
  269.   DIV32IA;
  270.   DIVA(25);
  271.    WHILE I>3 DO
  272.     BEGIN
  273.       I:=I-2;
  274.       DIV32IP;
  275.       SUBA;
  276.       DIVA(25);
  277.     END;
  278.    SUB32A;
  279.    I:=I2;                { COMPUTE TERM 2 TO I2 PLACES }
  280.    DIV4IA;
  281.    DIVA(239);
  282.    DIVA(239);
  283.     WHILE I>3 DO
  284.      BEGIN
  285.        I:=I-2;
  286.        DIV4IP;
  287.        SUBA;
  288.        DIVA(239);
  289.        DIVA(239);
  290.      END;
  291.        SUB4A;
  292.        DIVA(239);
  293.        SUBT;           { SUBTRACT TERM 2 FROM TERM 1 AND STORE IN A }
  294.        T2:=TIME;        { SET END OF COMPUTATION TIME }
  295. END;
  296.  
  297. PROCEDURE STORE;       { SAVE PI AND TIME TO  DISK FILES }
  298. VAR BUFF  : FILE OF INTEGER;
  299.         F : TEXT;
  300. BEGIN
  301. ASSIGN(BUFF,'PI.DTA');
  302. ASSIGN(F,'TIME.DTA');
  303. REWRITE(BUFF);
  304. REWRITE(F);
  305. FOR J:=0 TO M-2 DO
  306. WRITE(BUFF,A[J]);
  307.  
  308. WRITE(F,TI);           {Note, does not produce an ascii readable}
  309. WRITE(F,T2);           {output. It is for use by PILIST - code  }
  310.                        {for PILIST is at the end after PI5}
  311. CLOSE(BUFF);
  312. CLOSE(F);
  313. END;
  314.  
  315. procedure compute_time;
  316. var
  317.   h1,h2,m1,m2,s1,s2,hun1,hun2,code,x: integer;
  318.   s,temp: string[11];
  319. begin
  320.   val(copy(ti,1,2),h1,code);
  321.   val(copy(t2,1,2),h2,code);
  322.   val(copy(ti,4,2),m1,code);
  323.   val(copy(t2,4,2),m2,code);
  324.   val(copy(ti,7,2),s1,code);
  325.   val(copy(t2,7,2),s2,code);
  326.   val(copy(ti,10,2),hun1,code);
  327.   val(copy(t2,10,2),hun2,code);
  328.   if hun2<hun1 then begin
  329.      hun2:= hun2+100; s2:= s2-1; end;
  330.   if s2<s1 then begin
  331.      s2:= s2+60; m2:= m2-1; end;
  332.   if m2<m1 then begin
  333.      m2:= m2+60; h2:= h2-1; end;
  334.   if h2<h1 then begin
  335.      h2:= h2+24; end;
  336.   str(h2-h1:2,temp);
  337.   s:= temp+':';
  338.   str(m2-m1:2,temp);
  339.   s:= s+temp+':';
  340.   str(s2-s1:2,temp);
  341.   s:= s+temp+'.';
  342.   str(hun2-hun1:2,temp);
  343.   s:= s+temp;
  344.   for x:= 1 to length(s) do
  345.     if (s[x]=' ')and (s[x+1]='0') then begin
  346.       s[x+1]:= ' '; s[x+2]:= ' ' end;
  347.   writeln(s)
  348. {  writeln(h2-h1:2,':',m2-m1:2,':',s2-s1:2,'.',hun2-hun1:2); }
  349. end;
  350.  
  351. PROCEDURE PRINTPI;    { PRINT THE  FORMATED DIGITS OF PI }
  352. BEGIN
  353.    WRITE('PI=3.');
  354.    FOR J:=1 TO N   DO
  355.    BEGIN
  356. WRITE(A[J]);
  357.    IF J MOD 5 = 0 THEN WRITE(' ');
  358.    IF J MOD 50 = 0 THEN WRITE('  ',J:4,'  PL          ');
  359.    END;
  360.    WRITELN;WRITELN;
  361.  
  362.    WRITELN('ENDING   TIME = ',T2);
  363.    WRITELN('STARTING TIME = ',TI);
  364.    write  ('TOTAL TIME    = '); compute_time;
  365.    WRITELN;
  366. END;
  367.  
  368.  
  369. PROCEDURE HEADER;
  370.    BEGIN
  371.     WRITELN('                         THE COMPUTATION OF ',#227);
  372.     WRITELN('                  Press Control-Break to exit program. ');
  373.     WRITELN('                  ------------------------------------ ');
  374.     WRITELN;
  375.    END;
  376.  
  377. PROCEDURE START;               { PROMPT FOR INPUT AND INITIALIZE }
  378.    BEGIN
  379.      WRITELN('Input number of decimal places (Number < 2150)');
  380.       READLN(N);
  381.       TI:=TIME;
  382.       M:=N+2;
  383.       I2:=2*(3*M DIV 4)-4*TRUNC(LN(M))+5;
  384.    END;
  385.  
  386.  
  387. {MAIN PROGRAM}
  388. LABEL 20;
  389. BEGIN
  390. CLRSCR;
  391. HEADER;
  392. 20:START;
  393. COMPUTE;
  394. PRINTPI;
  395. {STORE;}                { REMOVE BRACKETS TO SAVE PI TO DISK }
  396. GOTO 20;                { DO AGAIN - PRESS BREAK TO EXIT PROGRAM }
  397. END.
  398.  
  399. { PROGRAM PILIST.PAS
  400. BLOCK WRITE TO PILIST.PAS AND DELETE FROM HERE WHEN YOU ARE SURE IT WORKS
  401. BE SURE TO TURN ON PRINTER.}
  402.  
  403. Uses Printer;
  404.  
  405. VAR
  406.    SIZE, J,K:INTEGER;
  407.    BUFF     :FILE OF INTEGER;
  408.    TBUFF    :TEXT;
  409.    T1, T2   :STRING[11];
  410. BEGIN
  411.       ASSIGN(BUFF,'PI.DTA');
  412.       ASSIGN(TBUFF,'TIME.DTA');
  413. RESET(BUFF);
  414. RESET(TBUFF);
  415. SIZE:=FILESIZE(BUFF)-1;
  416. WRITE('                 THE COMPUTATION OF PI TO ',SIZE,' PLACES');
  417. WRITELN;WRITELN;
  418. WRITE('PI=3.');
  419. WRITE(LST,'                 THE COMPUTATION OF PI TO ',SIZE,' PLACES');
  420. WRITELN(LST);WRITELN(LST);
  421. WRITE(LST,'PI=3.');
  422. READ(BUFF,J);
  423. K:=1;
  424. WHILE  NOT  EOF(BUFF)  DO
  425. BEGIN
  426. READ(BUFF,J);
  427. WRITE(J);
  428. WRITE(LST,J);
  429.    IF K MOD 5 = 0 THEN
  430.     BEGIN
  431.     WRITE(' ');
  432.     WRITE(LST,' ');
  433.     END;
  434.    IF K MOD 50 = 0 THEN
  435.    BEGIN
  436.     WRITE('  ',K:4,'  PL          ');
  437.     WRITE(LST,'  ',K:4,'  PL   ');
  438.     WRITELN(LST);
  439.     WRITE(LST,'     ');
  440.    END;
  441. K:=K+1;
  442. END;
  443.    WRITELN(LST);WRITELN(LST);
  444.    WRITELN;WRITELN;
  445.  
  446. READ(TBUFF,T1);
  447. READ(TBUFF,T2);
  448.  
  449.    WRITELN('ENDING   TIME = ',T2);
  450.    WRITELN('STARTING TIME = ',T1);
  451.    WRITELN(LST,'ENDING   TIME = ',T2);
  452.    WRITELN(LST,'STARTING TIME = ',T1);
  453.  
  454. END.
  455.